clear;
clc;

p=[0.52 0.50
    0.55 0.23
    0.50 0.27
    -0.46 -0.83
    0.51 0.80
    0.58 0.47
    -0.09 -0.26
    1.47 0.94
    -0.45 -0.21
    0.55 0.44
    0.23 0.05
    0.56 0.53
    0.23 -0.53
    1.68 1.38
    0 -0.52
    0.5 0.14
    0.16 0.31
    0.15 -0.63
    0.19 -0.63
    0.17 -0.75 
    1.22 0.67
    1.53 0.96
    1.53 1.01
    0.77 0.83
    0.47 0.46
    0.54 0.52];


wparam = p(:,1);
bparam = p(:,2);

global theta e12 eg ey

Nsim = 100000;
rng(10010);
ey = normrnd(0,1,Nsim,1);

param = wparam;
w.gamma1=param(1);  w.gamma2=param(2);
w.betay=param(3);   w.cy=param(4);      w.sigtheta=param(5);    
w.c1=param(6);      w.c1d=param(7);     w.phi1=param(8);    w.phi1d=param(9);    w.beta1=param(10);  w.beta1d=param(11);
w.c2=param(12);     w.c2d=param(13);    w.phi2=param(14);   w.phi2d=param(15);   w.beta2=param(16);  w.beta2d=param(17);
w.cg=param(18);     w.cr=param(19);     w.cr=param(20);
w.phig=param(21);   w.phir=param(22);   w.phim=param(23);
w.sigg=param(24);   w.sigm=param(25);   w.sigr=param(26);

param = bparam;
b.gamma1=param(1);  b.gamma2=param(2);
b.betay=param(3);   b.cy=param(4);      b.sigtheta=param(5);    
b.c1=param(6);      b.c1d=param(7);     b.phi1=param(8);    b.phi1d=param(9);    b.beta1=param(10);  b.beta1d=param(11);
b.c2=param(12);     b.c2d=param(13);    b.phi2=param(14);   b.phi2d=param(15);   b.beta2=param(16);  b.beta2d=param(17);
b.cg=param(18);     b.cr=param(19);     b.cr=param(20);
b.phig=param(21);   b.phir=param(22);   b.phim=param(23);
b.sigg=param(24);   b.sigm=param(25);   b.sigr=param(26);


%% Distribution of \theta
b.e12 = mvnrnd([0,0],[1,0;0,1],Nsim);
eg =normrnd(0,1,Nsim,1);
b.eg = eg*b.sigg;
w.e12 = mvnrnd([0,0],[1,0;0,1],Nsim);
w.eg = eg*w.sigg;

% p10=norminv(0.1); p20 = norminv(0.2);p40=norminv(0.4);p50=norminv(0.5);p60=norminv(0.6);p80=norminv(0.8);p90=norminv(0.9);
% 10th percentile
% thetao = repmat([p10,p20,p40,p50,p60,p80,p90],Nsim,1);
thetao = normrnd(0,1,Nsim,1);
% w.theta=thetao*w.sigtheta;
% b.theta=thetao*b.sigtheta;
btheta = b.cy+thetao*b.sigtheta;
wtheta = w.cy+thetao*w.sigtheta;


%% figure c4-1

e12=b.e12;
eg= b.eg;
theta = thetao*b.sigtheta;

% black teacher - black student
a=sim_all(bparam,0);
bias1_black=a(:,4);
bias2_black=a(:,5);

% white teacher - black student
a=sim_all(bparam,1);
bias1_nblack=a(:,4);
bias2_nblack=a(:,5);

e12=w.e12;
eg= w.eg;
theta = thetao*w.sigtheta;

% white teacher - white student
a=sim_all(wparam,0);
bias1_white=a(:,4);
bias2_white=a(:,5);

% white teacher - black student
a=sim_all(wparam,1);
bias1_nwhite=a(:,4);
bias2_nwhite=a(:,5);

figure
subplot(2,2,1);
data=[bias2_black btheta];
  [bandwidth,density,X,Y]=kde2d(data);

contour3(Y,X,density,100)
view([0 90]);
ylim([-1.1,1.1]);
  xlim([-3,3]);
  ylabel('Bias')
  xlabel('c_y+\theta')
  zlabel('density')
  title('Math Teacher Bias, Blacks')

subplot(2,2,2);
data=[bias2_nblack btheta];
  [bandwidth,density,X,Y]=kde2d(data);

contour3(Y,X,density,100)
view([0 90]);
ylim([-1.1,1.1]);
  xlim([-3,3]);
  ylabel('Bias')
  xlabel('c_y+\theta')
  zlabel('density')
  title('Math Teacher Bias, Other Race, Blacks')

  
subplot(2,2,3);
data=[bias2_white wtheta];
  [bandwidth,density,X,Y]=kde2d(data);

contour3(Y,X,density,100)
view([0 90]);
ylim([-1.1,1.1]);
  xlim([-3,3]);
  ylabel('Bias')
  xlabel('c_y+\theta')
  zlabel('density')
  title('Math Teacher Bias, Whites')

subplot(2,2,4);
data=[bias2_nwhite wtheta];
  [bandwidth,density,X,Y]=kde2d(data);

contour3(Y,X,density,100)
view([0 90]);
ylim([-1.1,1.1]);
  xlim([-3,3]);
  ylabel('Bias')
  xlabel('c_y+\theta')
  zlabel('density')
  title('Math Teacher Bias, Other Race, Whites')
saveas(gcf,'contour.jpg')

%  figure 2
[f_white_bias1_s,x_white_bias1_s]=ksdensity(bias1_white);
[f_white_bias2_s,x_white_bias2_s]=ksdensity(bias2_white);

[f_black_bias1_s,x_black_bias1_s]=ksdensity(bias1_black);
[f_black_bias2_s,x_black_bias2_s]=ksdensity(bias2_black);

[f_white_bias1_o,x_white_bias1_o]=ksdensity(bias1_nwhite);
[f_white_bias2_o,x_white_bias2_o]=ksdensity(bias2_nwhite);

[f_black_bias1_o,x_black_bias1_o]=ksdensity(bias1_nblack);
[f_black_bias2_o,x_black_bias2_o]=ksdensity(bias2_nblack);



figure
subplot(2,1,1);
plot(x_white_bias1_s,f_white_bias1_s,'k-',x_black_bias1_s,f_black_bias1_s,'k.-');
hold on;
plot([mean(bias1_white),mean(bias1_white)],ylim,'k--',[mean(bias1_black),mean(bias1_black)],ylim,'k:');
legend('Whites','Blacks');
xlabel('Bias');
ylabel('density');
title('Same Race ELA Teacher');

subplot(2,1,2);
plot(x_white_bias1_o,f_white_bias1_o,'k-',x_black_bias1_o,f_black_bias1_o,'k.-');
hold on;
plot([mean(bias1_nwhite),mean(bias1_nwhite)],ylim,'k--',[mean(bias1_nblack),mean(bias1_nblack)],ylim,'k:');
legend('Whites','Blacks');
xlabel('Bias');
ylabel('density');
title('Other Race ELA Teacher');
saveas(gcf,'distbias_ELA.jpg')



figure
subplot(2,1,1);
plot(x_white_bias2_s,f_white_bias2_s,'k-',x_black_bias2_s,f_black_bias2_s,'k.-');
hold on;
plot([mean(bias2_white),mean(bias2_white)],ylim,'k--',[mean(bias2_black),mean(bias2_black)],ylim,'k:');
legend('Whites','Blacks');
xlabel('Bias');
ylabel('density');
title('Same Race Math Teacher');

subplot(2,1,2);
plot(x_white_bias2_o,f_white_bias2_o,'k-',x_black_bias2_o,f_black_bias2_o,'k.-');
hold on;
plot([mean(bias2_nwhite),mean(bias2_nwhite)],ylim,'k--',[mean(bias2_nblack),mean(bias2_nblack)],ylim,'k:');
legend('Whites','Blacks');
xlabel('Bias');
ylabel('density');
title('Other Race Math Teacher');
saveas(gcf,'distbias_Math.jpg')

% cdf of Pr(Y=1), blacks, whites, same theta, same bias

theta = thetao*w.sigtheta;
a=sim_all(wparam,0);
[f_w,x_w]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);

theta = thetao*b.sigtheta;
a=sim_all(bparam,1);
[f_b,x_b]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);

% blacks, white \theta and G
theta = thetao*w.sigtheta;
bparam(4)=wparam(4);
bparam(21)=wparam(21);  bparam(24)=wparam(24);  bparam(18)=wparam(18);
a=sim_all(bparam,0);
[f_btg,x_btg]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);

% blacks, white \theta, G, and Ts
bparam(6)=wparam(6);bparam(8)=wparam(8);bparam(10)=wparam(10);
bparam(12)=wparam(12);bparam(14)=wparam(14);bparam(16)=wparam(16);
a=sim_all(bparam,0);
[f_btgt,x_btgt]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);



figure
plot(x_b,f_b,x_w,f_w,'k.-',x_btg,f_btg,x_btgt,f_btgt)
xlim([0,1]);ylim([0,1]);

% plot(x_b,f_b,'k',x_bdb,f_bdb,'k.',x_bdtsg,f_bdtsg,'k--',x_bdtdg,f_bdtdg,'k.-',x_w,f_w,'k-.')
%  legend('Blacks','Same \theta', 'Same Bias, diff GPA','Same Bias, Same GPA','Whites','location','southeast')

%% Pr(T)
wparam = p(:,1);
bparam = p(:,2);

theta_t = -1:0.01:1;

bg = b.cg+theta_t;
wg = w.cg+theta_t;

PrT1_1_bt_bs = normcdf(b.c1+b.phi1*theta_t+b.beta1*bg);
PrT1_1_bt_ws = normcdf(b.c1+b.c1d+(b.phi1+b.phi1d)*theta_t+(b.beta1+b.beta1d)*wg);
PrT1_1_btoe_bs = normcdf(b.c1+b.c1d+(b.phi1+b.phi1d)*theta_t+(b.beta1+b.beta1d)*bg);
PrT1_1_btwe_bs = normcdf(w.c1+w.phi1*theta_t+w.beta1*bg);
PrT1_1_wt_ws = normcdf(w.c1+w.phi1*theta_t+w.beta1*wg);
PrT1_1_wt_bs = normcdf(w.c1+w.c1d+(w.phi1+w.phi1d)*theta_t+(w.beta1+w.beta1d)*bg);
PrT1_1_wtoe_ws = normcdf(b.c1+b.c1d+(b.phi1+b.phi1d)*theta_t+(b.beta1+b.beta1d)*wg);
PrT1_1_wtbe_ws = normcdf(b.c1+b.phi1*theta_t+b.beta1*wg);

PrT2_1_bt_bs = normcdf(b.c2+b.phi2*theta_t+b.beta2*bg);
PrT2_1_bt_ws = normcdf(b.c2+b.c2d+(b.phi2+b.phi2d)*theta_t+(b.beta2+b.beta2d)*bg);
PrT2_1_btwe_bs = normcdf(w.c2+w.phi2*theta_t+w.beta2*bg);
PrT2_1_wt_ws = normcdf(w.c2+w.phi2*theta_t+w.beta2*wg);
PrT2_1_wt_bs = normcdf(w.c2+w.c2d+(w.phi2+w.phi2d)*theta_t+(w.beta2+w.beta2d)*wg);
PrT2_1_wtbe_ws = normcdf(b.c2+b.phi2*theta_t+b.beta2*wg);


plot(theta_t, PrT1_1_bt_bs, theta_t, PrT1_1_bt_ws, theta_t, PrT1_1_wt_ws, theta_t,PrT1_1_wt_bs, theta_t,PrT1_1_btwe_bs, theta_t,PrT1_1_wtbe_ws)
plot(theta_t, PrT2_1_bt_bs, theta_t, PrT2_1_bt_ws, theta_t, PrT2_1_wt_ws, theta_t,PrT2_1_wt_bs, theta_t,PrT2_1_btwe_bs, theta_t,PrT2_1_wtbe_ws)


plot(theta_t, PrT1_1_bt_bs, theta_t, PrT1_1_bt_ws)



